################################
## Plots for Muslim Attitudes ##
################################

######### Set working directory 

######### Load packages
library(readstata13)
library(ggplot2)
library(ggthemes)
library(grid)
library(gridExtra)
library(psych)
library(car)
library(Hmisc)
library(dotwhisker)
library(dplyr)
library(tidyr)
library(broom)
library(egg)
library(factoextra)

#########################
## Setting up the Data ##
#########################

data<-read.csv("Williamson_JEPS_Figure1_Data.csv")

head(data)

######### Multiplot function
# This function permits multiple plots to share the same legend. 

grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x + theme(legend.position="none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl),
                                            legend,
                                            ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend,
                                           ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  
  grid.newpage()
  grid.draw(combined)
  
  # return gtable invisibly
  invisible(combined)
  
}

scaleFUN <- function(x) sprintf("%.1f", x)

#######################
## Plot for Figure 1 ##
#######################

######### Main Effect Plots

# PCA Component Outcome
pca <- ggplot(data, aes(information, pca_attitudes, fill = as.factor(information))) + 
  stat_summary(fun.y = mean, geom = "bar",width=.75) + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 2),
               width = .4, alpha = 1,size=1.5) +
  theme_bw() + scale_fill_manual(values = c("0" = "grey85", "1" = "grey55"),
                                 labels = c("Control", "Treatment"),
                                 name = c("")) + 
  scale_y_continuous(labels=scaleFUN) +
  theme(axis.text.x = element_blank(), legend.text = element_text(size = 25),
        axis.ticks.x = element_blank(), legend.position = "bottom",
        axis.text.y = element_text(size = 20)) + xlab("") + ylab("") +
  ggtitle("PCA Component") + theme(text = element_text(size=18,face="bold")) + coord_cartesian(ylim=c(-.3, .3))
pca

# Thermometer Outcome
therm <- ggplot(data, aes(information, therm, fill = as.factor(information))) + 
  stat_summary(fun.y = mean, geom = "bar",width=.75) + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 2),
               width = .4, alpha = 1,size=1.5) +
  theme_bw() + scale_fill_manual(values = c("0" = "grey85", "1" = "grey55"),
                                 labels = c("Control", "Treatment"),
                                 name = c("")) + 
  scale_y_continuous(labels=scaleFUN) +
  theme(axis.text.x = element_blank(), legend.text = element_text(size = 25),
        axis.ticks.x = element_blank(), legend.position = "bottom",
        axis.text.y = element_text(size = 20)) + xlab("") + ylab("") +
  ggtitle("Thermometer") + theme(text = element_text(size=18,face="bold")) + coord_cartesian(ylim=c(50, 70))
therm

# Patriotism Outcome
pat <- ggplot(data, aes(information, patriotism, fill = as.factor(information))) + 
  stat_summary(fun.y = mean, geom = "bar",width=.75) + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 2),
               width = .4, alpha = 1,size=1.5) +
  theme_bw() + scale_fill_manual(values = c("0" = "grey85", "1" = "grey55"),
                                 labels = c("Control", "Treatment"),
                                 name = c("")) + 
  scale_y_continuous(labels=scaleFUN) +
  theme(axis.text.x = element_blank(), legend.text = element_text(size = 25),
        axis.ticks.x = element_blank(), legend.position = "bottom",
        axis.text.y = element_text(size = 20)) + xlab("") + ylab("") +
  ggtitle("Patriotism") + theme(text = element_text(size=18,face="bold")) + coord_cartesian(ylim=c(.5, .75))
pat

# Surveillance Outcome
surv <- ggplot(data, aes(information, surveillance, fill = as.factor(information))) + 
  stat_summary(fun.y = mean, geom = "bar",width=.75) + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 2),
               width = .4, alpha = 1,size=1.5) +
  theme_bw() + scale_fill_manual(values = c("0" = "grey85", "1" = "grey55"),
                                 labels = c("Control", "Treatment"),
                                 name = c("")) + 
  scale_y_continuous(labels=scaleFUN) +
  theme(axis.text.x = element_blank(), legend.text = element_text(size = 25),
        axis.ticks.x = element_blank(), legend.position = "bottom",
        axis.text.y = element_text(size = 20)) + xlab("") + ylab("") +
  ggtitle("Surveillance") + theme(text = element_text(size=18,face="bold")) + coord_cartesian(ylim=c(2, 3))
surv

# Muslim Refugee Ban Outcome
ban <- ggplot(data, aes(information, refugee_ban, fill = as.factor(information))) + 
  stat_summary(fun.y = mean, geom = "bar",width=.75) + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 2),
               width = .4, alpha = 1,size=1.5) +
  theme_bw() + scale_fill_manual(values = c("0" = "grey85", "1" = "grey55"),
                                 labels = c("Control", "Treatment"),
                                 name = c("")) + 
  scale_y_continuous(labels=scaleFUN) +
  theme(axis.text.x = element_blank(), legend.text = element_text(size = 25),
        axis.ticks.x = element_blank(), legend.position = "bottom",
        axis.text.y = element_text(size = 20)) + xlab("") + ylab("") +
  ggtitle("Refugee Ban") + theme(text = element_text(size=18,face="bold")) + coord_cartesian(ylim=c(2, 3))
ban

# Registration Outcome
reg <- ggplot(data, aes(information, registration, fill = as.factor(information))) + 
  stat_summary(fun.y = mean, geom = "bar",width=.75) + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 2),
               width = .4, alpha = 1,size=1.5) +
  theme_bw() + scale_fill_manual(values = c("0" = "grey85", "1" = "grey55"),
                                 labels = c("Control", "Treatment"),
                                 name = c("")) + 
  scale_y_continuous(labels=scaleFUN) +
  theme(axis.text.x = element_blank(), legend.text = element_text(size = 25),
        axis.ticks.x = element_blank(), legend.position = "bottom",
        axis.text.y = element_text(size = 20)) + xlab("") + ylab("") +
  ggtitle("Registration") + theme(text = element_text(size=18,face="bold")) + coord_cartesian(ylim=c(2, 3))
reg

# Combined Plots for Figure 1
grid_arrange_shared_legend(pca, therm, pat, surv, ban, reg, ncol =3, nrow = 2)

#######################
## Plot for Figure 2 ##
#######################

######### Robustness Check Plots

# Data for terrorism prime plots
hetero_terror<-read.csv("robust_data_terror.csv")
head(hetero_terror)
# Reordering factor variables
hetero_terror$model<-factor(hetero_terror$model,levels=c("No Prime", "With Prime", "Difference"))
hetero_terror$term<-factor(hetero_terror$term,levels=c("PCA","Thermometer","Patriotism","Surveillance","Muslim Ban","Registration"))
# Terror Prime Plot
terror_plot<-small_multiple(hetero_terror, dot_args = list(size=.75,color="black"),whisker_args = list(size = .1,color="black")) +
  theme_bw() +  geom_hline(yintercept = 0, colour = "grey60", linetype = 2,size=.75) +
  ggtitle("Priming Terrorism Fears") + 
  theme(plot.title = element_text(size=12,face="bold"),legend.position="none",
        axis.text.x=element_text(size=12,face="bold",angle = 60, hjust = 1),
        axis.text.y=element_text(size=10,face="bold"),
        strip.text.y=element_text(size=11,face="bold"),
        strip.background = element_rect(fill="grey90"))

# Data for non-PC prime plots:
hetero_pc<-read.csv("robust_data_sd.csv")
head(hetero_pc)
# Reordering factor variables
hetero_pc$model<-factor(hetero_pc$model,levels=c("No License", "With License", "Difference"))
hetero_pc$term<-factor(hetero_pc$term,levels=c("PCA","Thermometer","Patriotism","Surveillance","Muslim Ban","Registration"))
# Non-PC License Plot
pc_plot<-small_multiple(hetero_pc, dot_args = list(size=.75,color="black"),whisker_args = list(size = .1,color="black")) +
  theme_bw() +  geom_hline(yintercept = 0, colour = "grey60", linetype = 2,size=.75) +
  ggtitle("Licensing Non-PC Views") + 
  theme(plot.title = element_text(size=12,face="bold"),legend.position="none",
        axis.text.x=element_text(size=12,face="bold",angle = 60, hjust = 1),
        axis.text.y=element_text(size=10,face="bold"),
        strip.text.y=element_text(size=11,face="bold"),
        strip.background = element_rect(fill="grey90"))

# Combining Plots for Figure 2
ggarrange(terror_plot,pc_plot,ncol=2,nrow=1)

#######################
## Plot for Figure 3 ##
#######################

######### Heterogeneous Effect Plots

# Data for Elderly Plot
hetero_elderly<-read.csv("hetero_coef_data_elderly.csv")
head(hetero_elderly)
# Reordering Factor Variables
hetero_elderly$model<-factor(hetero_elderly$model,levels=c("Interaction","Elderly","Treatment"))
hetero_elderly$term<-factor(hetero_elderly$term,levels=c("PCA","Thermometer","Patriotism","Surveillance","Muslim Ban","Registration"))
# Elderly Plot 
elderly_plot<-small_multiple(hetero_elderly, dot_args = list(size=.75,color="black"),whisker_args = list(size = .1,color="black")) +
  theme_bw() +  geom_hline(yintercept = 0, colour = "grey60", linetype = 2,size=.75) +
  ggtitle("Elderly") + 
  theme(plot.title = element_text(size=12,face="bold"),legend.position="none",
        axis.text.x=element_text(size=12,face="bold",angle = 60, hjust = 1),
        axis.text.y=element_text(size=10,face="bold"),
        strip.text.y=element_text(size=11,face="bold"),
        strip.background = element_rect(fill="grey90"))

# Data for White Plot
hetero_white<-read.csv("hetero_coef_data_white.csv")
head(hetero_white)
# Reordering Factor Variables
hetero_white$model<-factor(hetero_white$model,levels=c("Interaction","White","Treatment"))
hetero_white$term<-factor(hetero_white$term,levels=c("PCA","Thermometer","Patriotism","Surveillance","Muslim Ban","Registration"))
# White Plot
white_plot<-small_multiple(hetero_white, dot_args = list(size=.75,color="black"),whisker_args = list(size = .1,color="black")) +
  theme_bw() +  geom_hline(yintercept = 0, colour = "grey60", linetype = 2,size=.75) +
  ggtitle("White") + 
  theme(plot.title = element_text(size=12,face="bold"),legend.position="none",
        axis.text.x=element_text(size=12,face="bold",angle = 60, hjust = 1),
        axis.text.y=element_text(size=10,face="bold"),
        strip.text.y=element_text(size=11,face="bold"),
        strip.background = element_rect(fill="grey90"))

# Data for No Muslim Plot
hetero_no<-read.csv("hetero_coef_data_no_muslim.csv")
head(hetero_no)
# Reordering Factor Variables
hetero_no$model<-factor(hetero_no$model,levels=c("Interaction","No Muslim","Treatment"))
hetero_no$term<-factor(hetero_no$term,levels=c("PCA","Thermometer","Patriotism","Surveillance","Muslim Ban","Registration"))
# No Muslim Plot
no_muslim_plot<-small_multiple(hetero_no, dot_args = list(size=.75,color="black"),whisker_args = list(size = .1,color="black")) +
  theme_bw() +  geom_hline(yintercept = 0, colour = "grey60", linetype = 2,size=.75) +
  ggtitle("No Muslim") + 
  theme(plot.title = element_text(size=12,face="bold"),legend.position="none",
        axis.text.x=element_text(size=12,face="bold",angle = 60, hjust = 1),
        axis.text.y=element_text(size=10,face="bold"),
        strip.text.y=element_text(size=11,face="bold"),
        strip.background = element_rect(fill="grey90"))

# Data for Republicans Plot
hetero_rep<-read.csv("hetero_coef_data_republican.csv")
head(hetero_rep)
# Reordering Factor Variables
hetero_rep$model<-factor(hetero_rep$model,levels=c("Interaction","Republican","Treatment"))
hetero_rep$term<-factor(hetero_rep$term,levels=c("PCA","Thermometer","Patriotism","Surveillance","Muslim Ban","Registration"))
# Republican Plot
rep_plot<-small_multiple(hetero_rep, dot_args = list(size=.75,color="black"),whisker_args = list(size = .1,color="black")) +
  theme_bw() +  geom_hline(yintercept = 0, colour = "grey60", linetype = 2,size=.75) +
  ggtitle("Republican") + 
  theme(plot.title = element_text(size=12,face="bold"),legend.position="none",
        axis.text.x=element_text(size=12,face="bold",angle = 60, hjust = 1),
        axis.text.y=element_text(size=10,face="bold"),
        strip.text.y=element_text(size=11,face="bold"),
        strip.background = element_rect(fill="grey90"))

# Combining Plots for Figure 3
ggarrange(elderly_plot,white_plot,no_muslim_plot,rep_plot,ncol=4,nrow=1)

############################
## Scree Plot in Appendix ##
############################

eigen<-read.csv("eigenvalues.csv")
head(eigen)

ggplot(data=eigen, aes(x=Component, y=Eigenvalue)) + theme_bw() +
  geom_line(lwd=1.5)+
  geom_point(size=5) + 
  geom_hline(yintercept = 1, colour = "grey60", linetype = 2,size=.75) +
  theme(axis.text.x=element_text(size=15,face="bold"),
        axis.text.y=element_text(size=15,face="bold"),
        axis.title=element_text(size=15,face="bold"))
